
      module mo_strato_rates
!=======================================================================
! ROUTINE
!   ratecon_sfstrat.f
!
!  Date...
!  15 August 2002
!  11 April  2008
!  15 December 2014
!
!  Programmed by...
!   Douglas E. Kinnison
!
! DESCRIPTION
!
! Derivation of the rate constant for reactions on
!   sulfate, NAT, and ICE aerosols.
!
!
! Sulfate Aerosol Reactions               Rxn#   Gamma
!   N2O5   + H2O(l)     =>  2HNO3         (1)    f(wt%)
!   ClONO2 + H2O(l)     =>  HOCl + HNO3   (2)    f(T,P,HCl,H2O,r)
!   BrONO2 + H2O(l)     =>  HOBr + HNO3   (3)    f(T,P,H2O,r)
!   ClONO2 + HCl(l)     =>  Cl2  + HNO3   (4)    f(T,P,HCl,H2O,r)
!   HOCl   + HCl(l)     =>  Cl2  + H2O    (5)    f(T,P,HCl,HOCl,H2O,r)
!   HOBr   + HCl(l)     =>  BrCl + H2O    (6)    f(T,P,HCl,HOBr,H2O,r)
!
! Nitric Acid Tri-hydrate Reactions       Rxn#    Gamma   Reference
!   N2O5   + H2O(s)     =>  2HNO3         (7)     4e-4   JPL10-6
!   ClONO2 + H2O(s)     =>  HOCl + HNO3   (8)     4e-3   JPL10-6
!   ClONO2 + HCl(s)     =>  Cl2  + HNO3   (9)     0.2    JPL10-6
!   HOCl   + HCl(s)     =>  Cl2  + H2O    (10)    0.1    JPL10-6
!   BrONO2 + H2O(s)     =>  HOBr + HNO3   (11)    0.006  Davies et JGR, 2003    
!
! WATER-ICE Aersol Reactions              Rxn#    Gamma
!   N2O5   + H2O(s)     =>  2HNO3         (12)     0.02   JPL10-6
!   ClONO2 + H2O(s)     =>  HOCl + HNO3   (13)     0.3    JPL10-6
!   BrONO2 + H2O(s)     =>  HOBr + HNO3   (14)     0.3    JPL10-6
!   ClONO2 + HCl(s)     =>  Cl2  + HNO3   (15)     0.3    JPL10-6
!   HOCl   + HCl(s)     =>  Cl2  + H2O    (16)     0.2    JPL10-6
!   HOBr   + HCl(s)     =>  BrCl + H2O    (17)     0.3    JPL10-6
!
! NOTE: The rate constants derived from species reacting with H2O are
!       first order (i.e., sec-1 units) - an example is N2O5 + H2O = 2HNO3.
!       Other reactions, e.g., ClONO2 + HCl have rate constants that
!       are second order (i.e., cm+3 molecules-1 sec-1 units). In all
!       of these types of reactions the derived first order rate constant
!       {0.25*(mean Velocity)*SAD*gamma} is divided by the HCl abundance
!       to derive the correct second order units.
!
! NOTE: Liquid Sulfate Aerosols...
!       See coding for references on how the Sulfate Aerosols were handled.
!       Approach follows Shi et al., JGR, 106, D20, 24259, 2001.
!
!
! INPUT:
!  ad      .    .... air density, molec. cm-3
!  pmid        ..... pressures, hPa
!  temp        ..... temperatures, K
!  rad_sulfate ..... Surface area density, cm2 cm-3
!  sad_sulfate ..... Surface area density, cm2 cm-3
!  sad_nat     ..... Surface area density, cm2 cm-3
!  sad_ice     ..... Surface area density, cm2 cm-3
!  brono2mv    ..... BrONO2 Volume Mixing Ratio
!  clono2mvr   ..... ClONO2 Volume Mixing Ratio
!  h2omvr      ..... H2O Volume Mixing Ratio
!  hclmvr      ..... HCl Volume Mixing Ratio
!  hobrmvr     ..... HOBr Volume Mixing Ratio
!  hoclmvr     ..... HOCl Volume Mixing Ratio
!  n2o5mvr     ..... N2O5 Volume Mixing Ratio
!
! OUTPUT:
!
!  rxt         ..... Rate constant (s-1 and cm3 sec-1 molec-1)
!=======================================================================

      private
      public :: ratecon_sfstrat, init_strato_rates, has_strato_chem

      integer :: id_brono2, id_clono2, id_hcl, id_hocl, &
           id_hobr, id_n2o5
      integer :: rid_het1,  rid_het2,  rid_het3,  rid_het4,  rid_het5, &
           rid_het6,  rid_het7,  rid_het8,  rid_het9,  rid_het10, &
           rid_het11, rid_het12, rid_het13, rid_het14, rid_het15, &
           rid_het16, rid_het17

      logical :: has_strato_chem 

      contains

        subroutine init_strato_rates

          use mo_chem_utls,      only : get_rxt_ndx, get_spc_ndx
          use mo_aero_settling,  only : strat_aer_settl_init
          use ppgrid,            only : pcols, pver

          implicit none

          integer :: ids(23)

          rid_het1  = get_rxt_ndx( 'het1' )
          rid_het2  = get_rxt_ndx( 'het2' )
          rid_het3  = get_rxt_ndx( 'het3' )
          rid_het4  = get_rxt_ndx( 'het4' )
          rid_het5  = get_rxt_ndx( 'het5' )
          rid_het6  = get_rxt_ndx( 'het6' )
          rid_het7  = get_rxt_ndx( 'het7' )
          rid_het8  = get_rxt_ndx( 'het8' )
          rid_het9  = get_rxt_ndx( 'het9' )
          rid_het10 = get_rxt_ndx( 'het10' )
          rid_het11 = get_rxt_ndx( 'het11' )
          rid_het12 = get_rxt_ndx( 'het12' )
          rid_het13 = get_rxt_ndx( 'het13' )
          rid_het14 = get_rxt_ndx( 'het14' )
          rid_het15 = get_rxt_ndx( 'het15' )
          rid_het16 = get_rxt_ndx( 'het16' )
          rid_het17 = get_rxt_ndx( 'het17' )

          id_brono2 = get_spc_ndx( 'BRONO2' )
          id_clono2 = get_spc_ndx( 'CLONO2' )
          id_hcl    = get_spc_ndx( 'HCL' )
          id_hocl   = get_spc_ndx( 'HOCL' )
          id_hobr   = get_spc_ndx( 'HOBR' )
          id_n2o5   = get_spc_ndx( 'N2O5' )

          ids(:) = (/ rid_het1, rid_het2, rid_het3, rid_het4, rid_het5, rid_het6, rid_het7, rid_het8, &
               rid_het9, rid_het10, rid_het11, rid_het12, rid_het13, rid_het14, rid_het15, &
               rid_het16, rid_het17, id_brono2, id_clono2, id_hcl, id_hocl, id_hobr, id_n2o5 /)

          has_strato_chem = all( ids(:) > 0 )

          if (.not. has_strato_chem) return

          call strat_aer_settl_init

        endsubroutine init_strato_rates

      subroutine ratecon_sfstrat( ncol, ad, pmid, temp, rad_sulfate, sad_sulfate, &
                                  sad_nat, sad_ice, h2ovmr, vmr, rxt, &
                                  gprob_n2o5, gprob_cnt_hcl, gprob_cnt_h2o, gprob_bnt_h2o, &
                                  gprob_hocl_hcl, gprob_hobr_hcl, wtper  )

      use shr_kind_mod, only : r8 => shr_kind_r8
      use chem_mods,    only : adv_mass, rxntot, gas_pcnst
      use ppgrid,       only : pcols, pver
      use mo_sad,       only : sad_top
      use cam_logfile,  only : iulog

      implicit none

!-----------------------------------------------------------------------
!	... dummy arguments
!-----------------------------------------------------------------------
      integer, intent(in) :: ncol                                 ! columns in chunk
      real(r8), dimension(ncol,pver,gas_pcnst), intent(in) :: &   ! species concentrations (mol/mol)
        vmr
      real(r8), dimension(ncol,pver), intent(in) :: &
        ad, &                                                     ! Air Density (molec. cm-3)
        rad_sulfate, &                                            ! Radius of Sulfate Aerosol (cm)
        sad_ice, &                                                ! ICE Surface Area Density (cm-1)
        sad_nat, &                                                ! NAT Surface Area Density (cm-1)
        sad_sulfate, &                                            ! Sulfate Surface Area Density (cm-1)
        h2ovmr                                                    ! water vapor volume mixing ratio( gas phase )
      real(r8), dimension(pcols,pver), intent(in) :: &
        pmid, &                                                   ! pressure (Pa)
        temp                                                      ! temperature (K)

      real(r8), intent(out) :: &
        rxt(ncol,pver,rxntot)                                     ! rate constants

      real(r8), dimension(ncol,pver), intent(out) :: &            ! diagnostics
        gprob_n2o5, &
        gprob_cnt_hcl, &
        gprob_cnt_h2o, &
        gprob_bnt_h2o, &
        gprob_hocl_hcl, &
        gprob_hobr_hcl, &
        wtper

!-----------------------------------------------------------------------
!  	... local variables
!-----------------------------------------------------------------------
        real(r8), parameter :: small_div  = 1.e-16_r8      ! for divid by excess species
	real(r8), parameter :: av_const   = 2.117265e4_r8  ! (8*8.31448*1000 / PI)
	real(r8), parameter :: pa2mb      = 1.e-2_r8       ! Pa to mb
	real(r8), parameter :: m2cm       = 100._r8        ! meters to cms

      integer :: &
        i, &                      ! altitude loop index
        k, &                      ! level loop index
        m                         ! species index

!-----------------------------------------------------------------------
!   	... variables for gamma calculations
!-----------------------------------------------------------------------
      real(r8) :: &
        brono2vmr, &                            ! BrONO2 Volume Mixing Ratio
        clono2vmr, &                            ! ClONO2 Volume Mixing Ratio
        hclvmr, &                               ! HCl Volume Mixing Ratio
        hcldeni, &                              ! inverse of HCl density
        cntdeni, &                              ! inverse of ClONO2 density
        hocldeni, &                             ! inverse of HOCl density
        hobrdeni, &                             ! inverse of HOBr density
        hoclvmr, &                              ! HOCl Volume Mixing Ratio
        hobrvmr, &                              ! HOBr Volume Mixing Ratio
        n2o5vmr                                 ! N2O5 Volume Mixing Ratio

        real(r8) :: &
        av_n2o5, &                              ! N2O5 Mean Velocity (cm s-1)
        av_clono2, &                            ! ClONO2 Mean Velocity (cm s-1)
        av_brono2, &                            ! BrONO2Mean Velocity (cm s-1)
        av_hocl, &                              ! HOCl Mean Velocity (cm s-1)
        av_hobr                                 ! HOBr Mean Velocity (cm s-1)

      real(r8) :: &
        pzero_h2o, &                            ! H2O sat vapor press (mbar)
        e0, e1, e2, e3, &                       ! coefficients for H2O sat vapor press.
        aw, &                                   ! Water activity
        m_h2so4, &                              ! H2SO4 molality (mol/kg)
        wt, &                                   ! wt % H2SO4
        y1, y2, &                               ! used in H2SO4 molality
        &  a1, b1, c1, d1, &                    ! used in H2SO4 molality
        a2, b2, c2, d2                          ! used in H2SO4 molality

        real(r8) :: &
        z1, z2, z3, &                           ! used in H2SO4 soln density
        den_h2so4, &                            ! H2SO4 soln density, g/cm3
        mol_h2so4, &                            ! Molality of H2SO4, mol / kg
        molar_h2so4, &                          ! Molarity of H2SO4, mol / l
        x_h2so4, &                              ! H2SO4 mole fraction
        aconst, tzero, &                        ! used in viscosity of H2SO4
        vis_h2so4, &                            ! H2SO4 viscosity
        ah, &                                   ! Acid activity, molarity units
        term1,term2,term3,term4, &              ! used in ah
        term5,term6,term7,term0, &
        T_limit, &                              ! temporary variable for temp (185-260K range)
        T_limiti, &                             ! 1./T_limit
        T_limitsq, &                            ! sqrt( T_limit )
        rad_sulf, &                             ! temporary variable for sulfate radius (cm)
        sadsulf, &                              ! temporary variable for sulfate radius (cm)
        sadice, &                               ! temporary variable for sulfate radius (cm)
        sadnat                                  ! temporary variable for sulfate radius (cm)

      real(r8) :: &
        C_cnt, S_cnt, &                         ! used in H_cnt
        H_cnt, &                                ! Henry's law coeff. for ClONO2
        H_hcl, &                                ! Henry's law coeff. for HCl
        D_cnt, &
        k_hydr, &
        k_h2o, &
        k_h, &
        k_hcl, &
        rdl_cnt, &
        f_cnt, &
        M_hcl, &
        atmos

      real(r8) :: &
        Gamma_b_h2o, &
        Gamma_cnt_rxn, &
        Gamma_b_hcl, &
        Gamma_s, &
        Fhcl, &
        Gamma_s_prime, &
        Gamma_b_hcl_prime, &
        Gamma_b, &
        gprob_rxn, &
        gprob_tot, &
        gprob_cnt
        
      real(r8) :: &
        D_hocl, &
        k_hocl_hcl, &
        C_hocl, &
        S_hocl, &
        H_hocl, &
        Gamma_hocl_rxn, &
        rdl_hocl, &
        f_hocl

      real(r8) :: &
        h1, h2, h3, &
        alpha

      real(r8) :: &
        C_hobr, &
        D_hobr, &
        aa, bb, cc, dd, &
        k_hobr_hcl, &
        k_dl, &
        k_wasch, &
        H_hobr, &
        rdl_hobr, &
        Gamma_hobr_rxn, &
        f_hobr

      real(r8) :: &
        pmb,&					! Pressure, mbar (hPa)
        pH2O_atm,&				! Partial press. H2O (atm)
        pH2O_hPa,&				! Partial press. H2O (hPa)
        pHCl_atm,&				! Partial press. HCl (atm)
        pCNT_atm                                ! Partial press. ClONO2 (atm)

!-----------------------------------------------------------------------
!     	... Used in pzero h2o calculation
!-----------------------------------------------------------------------
      real(r8), parameter :: wt_e0 = 18.452406985_r8
      real(r8), parameter :: wt_e1 = 3505.1578807_r8
      real(r8), parameter :: wt_e2 = 330918.55082_r8
      real(r8), parameter :: wt_e3 = 12725068.262_r8

      real(r8) :: &
        wrk, tmp

      

      if (.not. has_strato_chem) return

!-----------------------------------------------------------------------
!     	... intialize rate constants
!-----------------------------------------------------------------------
      do k = 1,pver
         rxt(:,k,rid_het1) = 0._r8
         rxt(:,k,rid_het2) = 0._r8
         rxt(:,k,rid_het3) = 0._r8
         rxt(:,k,rid_het4) = 0._r8
         rxt(:,k,rid_het5) = 0._r8
         rxt(:,k,rid_het6) = 0._r8
         rxt(:,k,rid_het7) = 0._r8
         rxt(:,k,rid_het8) = 0._r8
         rxt(:,k,rid_het9) = 0._r8
         rxt(:,k,rid_het10) = 0._r8
         rxt(:,k,rid_het11) = 0._r8
         rxt(:,k,rid_het12) = 0._r8
         rxt(:,k,rid_het13) = 0._r8
         rxt(:,k,rid_het14) = 0._r8
         rxt(:,k,rid_het15) = 0._r8
         rxt(:,k,rid_het16) = 0._r8
         rxt(:,k,rid_het17) = 0._r8

         gprob_n2o5(:,k)    = 0._r8
         gprob_cnt_h2o(:,k) = 0._r8
         gprob_cnt_hcl(:,k) = 0._r8
         gprob_bnt_h2o(:,k) = 0._r8
         gprob_hocl_hcl(:,k)= 0._r8
         gprob_hobr_hcl(:,k)= 0._r8
         wtper(:,k)         = 0._r8
      end do

!-----------------------------------------------------------------------
!     	... set rate constants
!-----------------------------------------------------------------------
Level_loop : &
      do k = sad_top+1,pver
column_loop : &
         do i = 1,ncol
!-----------------------------------------------------------------------
!	... set species, pmb, and atmos
!-----------------------------------------------------------------------
	    brono2vmr    = vmr(i,k,id_brono2)
	    clono2vmr    = vmr(i,k,id_clono2)
	    hclvmr       = vmr(i,k,id_hcl)
	    hoclvmr      = vmr(i,k,id_hocl)
	    hobrvmr      = vmr(i,k,id_hobr)
	    if( hclvmr > 0._r8 ) then
	       hcldeni  = 1._r8/(hclvmr*ad(i,k))
	    end if
	    if( clono2vmr > 0._r8 ) then
	       cntdeni  = 1._r8/(clono2vmr*ad(i,k))
	    end if
	    if( hoclvmr > 0._r8 ) then
	       hocldeni  = 1._r8/(hoclvmr*ad(i,k))
	    end if
	    if( hobrvmr > 0._r8 ) then
	       hobrdeni  = 1._r8/(hobrvmr*ad(i,k))
	    end if
	    n2o5vmr      = vmr(i,k,id_n2o5)
            sadsulf      = sad_sulfate(i,k)
            sadnat       = sad_nat(i,k)
	    sadice       = sad_ice(i,k)
            pmb          = pa2mb*pmid(i,k)
            atmos        = pmb/1013.25_r8

!-----------------------------------------------------------------------
!  	... setup for stratospheric aerosols
!           data range set: 185K - 240K;    Tabazedeh GRL, 24, 1931, 1997
!-----------------------------------------------------------------------
            T_limit   = max( temp(i,k),185._r8 )
            T_limit   = min( T_limit,240._r8 )
            T_limiti  = 1._r8/T_limit
            T_limitsq = sqrt( T_limit )

!-----------------------------------------------------------------------
!     .... Average velocity (8RT*1000/(PI*MW))**1/2 * 100.(units cm s-1)
!     .... or (av_const*T/M2)**1/2
!-----------------------------------------------------------------------
	    wrk       = av_const*T_limit
            av_n2o5   = sqrt( wrk/adv_mass(id_n2o5) )*m2cm
            av_clono2 = sqrt( wrk/adv_mass(id_clono2) )*m2cm
            av_brono2 = sqrt( wrk/adv_mass(id_brono2) )*m2cm
            av_hocl   = sqrt( wrk/adv_mass(id_hocl) )*m2cm
            av_hobr   = sqrt( wrk/adv_mass(id_hobr) )*m2cm
has_sadsulf : &
            if( sadsulf > 0._r8 ) then
!-----------------------------------------------------------------------
!     .... Partial Pressure of H2O, ClONO2, and HCl in atmospheres
!-----------------------------------------------------------------------
               if( hclvmr > 0._r8 ) then
                  pHCl_atm  = hclvmr*atmos
               else
                  pHCl_atm  = 0._r8
               end if

               if( clono2vmr > 0._r8 ) then
                  pCNT_atm  = clono2vmr*atmos
               else
                  pCNT_atm  = 0._r8
               end if

               if( h2ovmr(i,k) > 0._r8 ) then
                  pH2O_atm  = h2ovmr(i,k)*atmos
               else
                  pH2O_atm  = 0._r8
               end if
!-----------------------------------------------------------------------
!     .... Partial Pressure of H2O in hPa
!-----------------------------------------------------------------------
               pH2O_hpa = h2ovmr(i,k)*pmb
!-----------------------------------------------------------------------
!     .... Calculate the h2so4 Wt% and Activity of H2O - 
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     ... Saturation Water Vapor Pressure (mbar)
!-----------------------------------------------------------------------
               pzero_h2o = exp( wt_e0 - T_limiti*(wt_e1 + T_limiti*(wt_e2 - T_limiti*wt_e3)) )
!-----------------------------------------------------------------------
!     ... H2O activity
!     ... if the activity of H2O goes above 1.0, wt% can go negative
!-----------------------------------------------------------------------
               aw = ph2o_hpa / pzero_h2o
               aw = min( aw,1._r8 )
               aw = max( aw,.001_r8 )
!-----------------------------------------------------------------------
!     ... h2so4 Molality (mol/kg)
!-----------------------------------------------------------------------
               if( aw <= .05_r8 ) then
                  a1 = 12.37208932_r8
                  b1 = -0.16125516114_r8
                  c1 = -30.490657554_r8
                  d1 = -2.1133114241_r8
                  a2 = 13.455394705_r8
                  b2 = -0.1921312255_r8
                  c2 = -34.285174607_r8
                  d2 = -1.7620073078_r8
               else if( aw > .05_r8 .and. aw < .85_r8 ) then
                  a1 = 11.820654354_r8
                  b1 = -0.20786404244_r8
                  c1 = -4.807306373_r8
                  d1 = -5.1727540348_r8
                  a2 = 12.891938068_r8
                  b2 = -0.23233847708_r8
                  c2 = -6.4261237757_r8
                  d2 = -4.9005471319_r8
               else
                  a1 = -180.06541028_r8
                  b1 = -0.38601102592_r8
                  c1 = -93.317846778_r8
                  d1 = 273.88132245_r8
                  a2 = -176.95814097_r8
                  b2 = -0.36257048154_r8
                  c2 = -90.469744201_r8
                  d2 = 267.45509988_r8
               end if
!-----------------------------------------------------------------------
!     ... h2so4 mole fraction
!-----------------------------------------------------------------------
               y1       = a1*(aw**b1) + c1*aw + d1
               y2       = a2*(aw**b2) + c2*aw + d2
               m_h2so4  = y1 + ((T_limit - 190._r8)*(y2 - y1)) / 70._r8
!-----------------------------------------------------------------------
!     ... h2so4 Weight Percent
!-----------------------------------------------------------------------
               wt = 9800._r8*m_h2so4 / (98._r8*m_h2so4  + 1000._r8)
               wtper(i,k) = wt
!-----------------------------------------------------------------------
!     .... Parameters for h2so4 Solution
!-----------------------------------------------------------------------
!     ... h2so4 Solution Density (g/cm3)
!-----------------------------------------------------------------------
	       wrk = T_limit*T_limit
               z1 =  .12364_r8  - 5.6e-7_r8*wrk
               z2 = -.02954_r8  + 1.814e-7_r8*wrk
               z3 =  2.343e-3_r8 - T_limit*1.487e-6_r8 - 1.324e-8_r8*wrk
!-----------------------------------------------------------------------
!     ... where mol_h2so4 is molality in mol/kg
!-----------------------------------------------------------------------
               den_h2so4 = 1._r8 + m_h2so4*(z1 + z2*sqrt(m_h2so4) + z3*m_h2so4)
!-----------------------------------------------------------------------
!     ... h2so4 Molarity, mol / l
!-----------------------------------------------------------------------
               molar_h2so4 = den_h2so4*wt/9.8_r8
!-----------------------------------------------------------------------
!     ... h2so4 Mole fraction
!-----------------------------------------------------------------------
               x_h2so4   = wt / (wt + ((100._r8 - wt)*98._r8/18._r8))
               term1     = .094_r8 - x_h2so4*(.61_r8 - 1.2_r8*x_h2so4)
               term2     = (8515._r8 - 10718._r8*(x_h2so4**.7_r8))*T_limiti
               H_hcl     = term1 * exp( -8.68_r8 + term2 )
               M_hcl     = H_hcl*pHCl_atm
!-----------------------------------------------------------------------
!     ... h2so4 solution viscosity
!-----------------------------------------------------------------------
               aconst    = 169.5_r8 + wt*(5.18_r8 - wt*(.0825_r8 - 3.27e-3_r8*wt))
               tzero     = 144.11_r8 + wt*(.166_r8 - wt*(.015_r8 - 2.18e-4_r8*wt))
               vis_h2so4 = aconst/(T_limit**1.43_r8) * exp( 448._r8/(T_limit - tzero) )
!-----------------------------------------------------------------------
!     ... Acid activity in molarity
!-----------------------------------------------------------------------
               term1 = 60.51_r8
               term2 = .095_r8*wt
	       wrk   = wt*wt
               term3 = .0077_r8*wrk
               term4 = 1.61e-5_r8*wt*wrk
               term5 = (1.76_r8 + 2.52e-4_r8*wrk) * T_limitsq
               term6 = -805.89_r8 + (253.05_r8*(wt**.076_r8))
               term7 = T_limitsq
               ah    = exp( term1 - term2 + term3 - term4 - term5 + term6/term7 )
	       if( ah <= 0._r8 ) then
	          write(iulog,*) 'ratecon: ah <= 0 at i,k, = ',i,k
	          write(iulog,*) 'ratecon: term1,term2,term3,term4,term5,term6,term7,wt,T_limit,ah = ', &
	                               term1,term2,term3,term4,term5,term6,term7,wt,T_limit,ah 
	       end if

	       wrk      = .25_r8*sadsulf
               rad_sulf = max( rad_sulfate(i,k),1.e-6_r8 )
!-----------------------------------------------------------------------
!     N2O5 + H2O(liq) =>  2.00*HNO3  Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
                  term0 = -25.5265_r8 - wt*(.133188_r8 - wt*(.00930846_r8 - 9.0194e-5_r8*wt))
                  term1 = 9283.76_r8 + wt*(115.345_r8 - wt*(5.19258_r8 - .0483464_r8*wt))
                  term2 = -851801._r8 - wt*(22191.2_r8 - wt*(766.916_r8 - 6.85427_r8*wt))
                  gprob_n2o5(i,k) = exp( term0 + T_limiti*(term1 + term2*T_limiti) )
                  rxt(i,k,rid_het1) = max( 0._r8,wrk*av_n2o5*gprob_n2o5(i,k) )

!-----------------------------------------------------------------------
!     ClONO2 + H2O(liq) =  HOCl + HNO3   Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
! 	... NOTE: Aerosol radius in units of cm.
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     	... Radius sulfate set (from sad module)
!           Set min radius to 0.01 microns (1e-6 cm)
!           Typical radius is 0.1 microns (1e-5 cm)
!           f_cnt may go negative under if not set.
!-----------------------------------------------------------------------
                  C_cnt         = 1474._r8*T_limitsq
                  S_cnt         = .306_r8 + 24._r8*T_limiti
                  term1         = exp( -S_cnt*molar_h2so4 )
                  H_cnt         = 1.6e-6_r8 * exp( 4710._r8*T_limiti )*term1
                  D_cnt         = 5.e-8_r8*T_limit / vis_h2so4
                  k_h           = 1.22e12_r8*exp( -6200._r8*T_limiti )
                  k_h2o         = 1.95e10_r8*exp( -2800._r8*T_limiti )
                  k_hydr        = (k_h2o + k_h*ah)*aw
                  k_hcl         = 7.9e11_r8*ah*D_cnt*M_hcl
                  rdl_cnt       = sqrt( D_cnt/(k_hydr + k_hcl) )
                  term1         = 1._r8/tanh( rad_sulf/rdl_cnt )
                  term2         = rdl_cnt/rad_sulf
                  f_cnt         = term1 - term2
                  if( f_cnt > 0._r8 ) then
                     term1         = 4._r8*H_cnt*.082_r8*T_limit
                     term2         = sqrt( D_cnt*k_hydr )
                     Gamma_b_h2o   = term1*term2/C_cnt
                     term1         = sqrt( 1._r8 + k_hcl/k_hydr )
                     Gamma_cnt_rxn = f_cnt*Gamma_b_h2o*term1
                     Gamma_b_hcl   = Gamma_cnt_rxn*k_hcl/(k_hcl + k_hydr)
                     term1         = exp( -1374._r8*T_limiti )
                     Gamma_s       = 66.12_r8*H_cnt*M_hcl*term1
		     if( pHCl_atm > 0._r8 ) then
                        term1      = .612_r8*(Gamma_s+Gamma_b_hcl)* pCNT_atm/pHCl_atm
                        Fhcl       = 1._r8/(1._r8 + term1)
		     else
                        Fhcl       = 1._r8
		     end if
                     Gamma_s_prime     = Fhcl*Gamma_s
                     Gamma_b_hcl_prime = Fhcl*Gamma_b_hcl
                     term1         = Gamma_cnt_rxn*k_hydr
                     term2         = k_hcl + k_hydr
                     Gamma_b       = Gamma_b_hcl_prime + (term1/term2)
                     term1         = 1._r8 / (Gamma_s_prime + Gamma_b)
                     gprob_cnt     = 1._r8 / (1._r8 + term1)
                     term1         = Gamma_s_prime + Gamma_b_hcl_prime
                     term2         = Gamma_s_prime + Gamma_b
                     gprob_cnt_hcl(i,k) = gprob_cnt * term1/term2
                     gprob_cnt_h2o(i,k) = gprob_cnt - gprob_cnt_hcl(i,k)
                  else
                     gprob_cnt_h2o(i,k) = 0._r8
                     gprob_cnt_hcl(i,k) = 0._r8
                     Fhcl          = 1._r8
                  end if

                  rxt(i,k,rid_het2) = max( 0._r8,wrk*av_clono2*gprob_cnt_h2o(i,k) )

!-----------------------------------------------------------------------
!  	... BrONO2 + H2O(liq) =  HOBr + HNO3   Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
                  h1    = 29.24_r8
                  h2    = -.396_r8
                  h3    = .114_r8
                  alpha = .805_r8
                  gprob_rxn = exp( h1 + h2*wt ) + h3
                  term1     = 1._r8/alpha
                  term2     = 1._r8/gprob_rxn
                  gprob_bnt_h2o(i,k) = 1._r8 / (term1 + term2)
                  rxt(i,k,rid_het3) = max( 0._r8,wrk*av_brono2*gprob_bnt_h2o(i,k) )

!-----------------------------------------------------------------------
!     	... ClONO2 + HCl(liq) =  Cl2  + HNO3  Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
               if( hclvmr > small_div .and. clono2vmr > small_div ) then
                 if ( hclvmr > clono2vmr ) then
                    rxt(i,k,rid_het4) = max( 0._r8,wrk*av_clono2*gprob_cnt_hcl(i,k) )*hcldeni
                 else
                    rxt(i,k,rid_het4) = max( 0._r8,wrk*av_clono2*gprob_cnt_hcl(i,k) )*cntdeni
                 end if
               end if

!-----------------------------------------------------------------------
!     	... HOCl + HCl(liq) =  Cl2 + H2O   Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     	... Radius sulfate set (from sad module)
!           Set min radius to 0.01 microns (1e-6 cm)
!           Typical radius is 0.1 microns (1e-5 cm)
!           f_hocl may go negative under if not set.
!-----------------------------------------------------------------------
	          if( pCNT_atm > 0._r8 ) then
                     D_hocl          = 6.4e-8_r8*T_limit/vis_h2so4
                     k_hocl_hcl      = 1.25e9_r8*ah*D_hocl*M_hcl
                     C_hocl          = 2009._r8*T_limitsq
                     S_hocl          = .0776_r8 + 59.18_r8*T_limiti
                     term1           = exp( -S_hocl*molar_h2so4 )
                     H_hocl          = 1.91e-6_r8 * exp( 5862.4_r8*T_limiti )*term1
                     term1           = 4._r8*H_hocl*.082_r8*T_limit
                     term2           = sqrt( D_hocl*k_hocl_hcl )
                     Gamma_hocl_rxn  = term1*term2/C_hocl
                     rdl_hocl        = sqrt( D_hocl/k_hocl_hcl )
                     term1           = 1._r8/tanh( rad_sulf/rdl_hocl )
                     term2           = rdl_hocl/rad_sulf
                     f_hocl          = term1 - term2
                     if( f_hocl > 0._r8 ) then
                        term1           = 1._r8 / (f_hocl*Gamma_hocl_rxn*Fhcl)
                        gprob_hocl_hcl(i,k)  = 1._r8 / (1._r8 + term1)
                     else
                        gprob_hocl_hcl(i,k)  = 0._r8
                     end if

                     if( hclvmr > small_div .and. hoclvmr > small_div ) then
                       if ( hclvmr > hoclvmr ) then
                         rxt(i,k,rid_het5) = max( 0._r8,wrk*av_hocl*gprob_hocl_hcl(i,k) )*hcldeni
                       else
                         rxt(i,k,rid_het5) = max( 0._r8,wrk*av_hocl*gprob_hocl_hcl(i,k) )*hocldeni
                       end if
                     end if
	          end if

!-----------------------------------------------------------------------
!     	... HOBr + HCl(liq) =  BrCl + H2O  Sulfate Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!   	... Radius sulfate set (from sad module)
!           Set min radius to 0.01 microns (1e-6 cm)
!           Typical radius is 0.1 microns (1e-5 cm)
!           f_hobr may go negative under if not set.
!-----------------------------------------------------------------------
                  C_hobr          = 1477._r8*T_limitsq
                  D_hobr          = 9.e-9_r8
!-----------------------------------------------------------------------
!     	...  Taken from Waschewsky and Abbat
!            Dave Hanson (PC) suggested we divide this rc by eight to agee
!            with his data (Hanson, 108, D8, 4239, JGR, 2003).
!            k1=k2*Mhcl for gamma(HOBr)
!-----------------------------------------------------------------------
                  k_wasch         = .125_r8 * exp( .542_r8*wt - 6440._r8*T_limiti + 10.3_r8)
!-----------------------------------------------------------------------
!     	... Taken from Hanson 2002.
!-----------------------------------------------------------------------
                  H_hobr          = exp( -9.86_r8 + 5427._r8*T_limiti )
                  k_dl            = 7.5e14_r8*D_hobr*2._r8                   ! or  7.5e14*D *(2nm)
!-----------------------------------------------------------------------
!  	... If k_wasch is GE than the diffusion limit...
!-----------------------------------------------------------------------
		  if( M_hcl > 0._r8 ) then
                     if( k_wasch >= k_dl ) then
                        k_hobr_hcl   = k_dl * M_hcl
                     else
                        k_hobr_hcl   = k_wasch * M_hcl
                     end if
                     term1           = 4._r8*H_hobr*.082_r8*T_limit
                     term2           = sqrt( D_hobr*k_hobr_hcl )
                     tmp             = rad_sulf/term2
                     Gamma_hobr_rxn  = term1*term2/C_hobr
                     rdl_hobr        = sqrt( D_hobr/k_hobr_hcl )
		     if( tmp < 1.e2_r8 ) then
                        term1           = 1._r8/tanh( rad_sulf/rdl_hobr )
		     else
                        term1           = 1._r8
		     end if
                     term2           = rdl_hobr/rad_sulf
                     f_hobr          = term1 - term2
                     if( f_hobr > 0._r8 ) then
                        term1              = 1._r8 / (f_hobr*Gamma_hobr_rxn)
                        gprob_hobr_hcl(i,k)= 1._r8 / (1._r8 + term1)
                     else
                        gprob_hobr_hcl(i,k)= 0._r8
                     end if
                     if( hclvmr > small_div .and. hobrvmr > small_div ) then
                       if ( hclvmr > hobrvmr ) then
                          rxt(i,k,rid_het6) = max( 0._r8,wrk*av_hobr*gprob_hobr_hcl(i,k) )*hcldeni
                       else
                          rxt(i,k,rid_het6) = max( 0._r8,wrk*av_hobr*gprob_hobr_hcl(i,k) )*hobrdeni    
                       end if           
                     end if
		  end if
 
            end if has_sadsulf

has_sadnat : &
	    if( sadnat > 0._r8 ) then
	       wrk = .25_r8*sadnat
!-----------------------------------------------------------------------
!     	... N2O5 + H2O(s) => 2HNO3  NAT Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     ... gprob based on JPL10-6 for NAT.
!         also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
!                 gprob_tot     = 4.e-4
!-----------------------------------------------------------------------
                rxt(i,k,rid_het7)  = wrk*av_n2o5*4.e-4_r8

!-----------------------------------------------------------------------
!     ClONO2 + H2O(s) => HNO3 + HOCl  NAT Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     ... gprob based on JPL10-6 for NAT.
!         also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
!                 gprob_tot    = 0.004
!-----------------------------------------------------------------------
                rxt(i,k,rid_het8) = wrk*av_clono2*4.0e-3_r8

!-----------------------------------------------------------------------
!     	... ClONO2 + HCl(s) => HNO3 + Cl2, NAT Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     ... gprob based on JPL10-6 for NAT.
!         also see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
!                 gprob_tot   = 0.2
!-----------------------------------------------------------------------
                if( hclvmr > small_div .and. clono2vmr > small_div ) then
                   if ( hclvmr > clono2vmr ) then
                      rxt(i,k,rid_het9) = wrk*av_clono2*0.2_r8*hcldeni
                   else
                      rxt(i,k,rid_het9) = wrk*av_clono2*0.2_r8*cntdeni  
                   end if
                end if

!-----------------------------------------------------------------------
!     	... HOCl + HCl(s) => H2O + Cl2  NAT Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     ... gprob based on JPL10-6 for NAT.
!         see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
!         and Abbatt and Molina, GRL, 19, 461-464, 1992.
!                 gprob_tot   = 0.1
!-----------------------------------------------------------------------
               if( hclvmr > small_div .and. hoclvmr > small_div ) then
                   if ( hclvmr > hoclvmr ) then
                      rxt(i,k,rid_het10) = wrk*av_hocl*0.1_r8*hcldeni
                   else
                      rxt(i,k,rid_het10) = wrk*av_hocl*0.1_r8*hocldeni
                   end if
               end if

!-----------------------------------------------------------------------
!     	... BrONO2 + H2O(s) => HOBr + HNO3  NAT Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!       ... Davies et al., 
!           JGR, 108, NO. D5, 8322, doi:10.1029/2001JD000445, 2003
!                 gprob_tot   = 0.006
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het11) = wrk*av_brono2*0.006_r8

            end if has_sadnat

has_sadice : &
	    if( sadice > 0._r8 ) then
	         wrk = .25_r8*sadice
!-----------------------------------------------------------------------
!     N2O5 + H2O(s) => 2HNO3  ICE Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!       ... gprob based on JPL10-6 for ICE.
!           also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
!                 gprob_tot    = 0.02
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het12) = wrk*av_n2o5*0.02_r8

!-----------------------------------------------------------------------
!     	... ClONO2 + H2O(s) => HNO3 + HOCl  ICE Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     	... gprob based on JPL10-6 for ICE.
!     	    also see Hanson and Ravi, JGR, 96, 17307-17314, 1991.
!                 gprob_tot    = 0.3
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het13) = wrk*av_clono2*0.3_r8

!-----------------------------------------------------------------------
!     	... BrONO2 + H2O(s) => HNO3 + HOBr  ICE Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!     	... gprob based on JPL10-6 for ICE.
!           also see Hanson and Ravi, JPC, 97, 2802-2803, 1993.
!           could be as high as 1.0
!                 gprob_tot    = 0.3
!-----------------------------------------------------------------------
                  rxt(i,k,rid_het14) = wrk*av_brono2*0.3_r8

!-----------------------------------------------------------------------
!     ClONO2 + HCl(s) => HNO3 + Cl2, ICE Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!       ... gprob based on JPL10-6 for ICE.
!           also see Hanson and Ravi, GRL, 15, 17-20, 1988.
!           also see Lue et al.,
!                 gprob_tot    = 0.3
!-----------------------------------------------------------------------
                 if( hclvmr > small_div .and. clono2vmr > small_div ) then
                     if ( hclvmr > clono2vmr ) then
                        rxt(i,k,rid_het15) = wrk*av_clono2*0.3_r8*hcldeni
                     else
                        rxt(i,k,rid_het15) = wrk*av_clono2*0.3_r8*cntdeni
                     end if
                 end if
!
!-----------------------------------------------------------------------
!     	... HOCl + HCl(s) => H2O + Cl2, ICE Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!       ... gprob based on JPL10-6 for ICE.
!           also see Hanson and Ravi, JPC, 96, 2682-2691, 1992.
!           also see Abbatt and Molina, GRL, 19, 461-464, 1992.
!                 gprob_tot   = 0.2
!-----------------------------------------------------------------------
                 if( hoclvmr > small_div .and. hclvmr > small_div ) then
                     if ( hclvmr > hoclvmr ) then
                        rxt(i,k,rid_het16) = wrk*av_hocl*0.2_r8*hcldeni
                     else
                        rxt(i,k,rid_het16) = wrk*av_hocl*0.2_r8*hocldeni
                     end if
                  end if

!-----------------------------------------------------------------------
!     HOBr + HCl(s) => H2O + BrCl, ICE Aerosol Reaction
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!       ... gprob based on JPL10-6 for ICE.
!           Abbatt GRL, 21, 665-668, 1994.
!                    gprob_tot   = 0.3
!-----------------------------------------------------------------------
                  if( hobrvmr > small_div .and. hclvmr > small_div ) then
                     if ( hclvmr > hobrvmr ) then
                        rxt(i,k,rid_het17) = wrk*av_hobr*0.3_r8*hcldeni
                     else
                        rxt(i,k,rid_het17) = wrk*av_hobr*0.3_r8*hobrdeni
                     end if
                  end if

            end if has_sadice
         end do column_loop
      end do Level_loop

      end subroutine ratecon_sfstrat

      end module mo_strato_rates
